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I .  INTRODUCTION 

The  millimeter  wave  radios  currently  under  development  for  use 
by  US  Army  personnel  may  have  to  transmit  information  through  groves 
of  trees  or  be  placed  within  a  wooded  area  for  camouflage  purposes. 
Recent  experimental  results  [VE,SJ]  have  indicated  that  such  communi¬ 
cation  capability  is  feasible.  To  support  this  conclusion  this 
study  was  initiated.  A  theory  of  millimeter-wave  propagation  in  veg¬ 
etation  (forests)  using  transport  theory  evolved. 

In  the  millimeter  range  a  forest  can  be  characterized  as  a 
random  distribution  of  particles  which  absorb  and  scatter  electro¬ 
magnetic  energy  in  all  directions.  Transport  theory  accounts  for 
multiple  scattering  effects  fully  but  neglects  interference  phenom¬ 
ena,  which  is  justified  by  the  experimental  results.  In  transport 
theory  the  total  field  intensity  is  separated  into  two  parts  -  the 
coherent  intensity  and  the  incoherent  intensity.  The  polarization 
properties  of  the  coherent  component  are  the  same  as  those  of  the 
incident  field,  which  is  taken  to  be  a  plane  wave  impinging  upon  the 
forest;  the  forest  is  assumed  to  be  a  slab  or  half-space  region. 

The  incoherent  intensity,  on  the  other  hand,  is  described  by  using 
the  four  Stokes'  parameters.  The  theory,  thus,  involves  a  system  of 
four  coupled  integro-dif ferential  equations  for  the  specific  inten¬ 
sity  of  the  incoherent  component. 

By  making  reasonable  assumption  concerning  the  statistical  scat¬ 
tering  characteristics  of  the  forest  medium,  it  was  found  that  for 
the  first  of  the  Stokes'  parameters  a  single,  uncoupled  equation  is 
obtained  which  is  identical  to  the  conventional  scalar  transport 
equation.  The  first  Stokes'  parameter  is  the  quantity  of  primary  in¬ 
terest  since  it  is,  by  definition,  the  total  intensity  of  radiation, 
i.e.,  it  represents  the  sum  of  the  intensities  of  two  orthogonal 
polarizations.  This  being  the  case,  only  the  scalar  transport  equa¬ 
tion  was  needed.  The  relevant  equations  are  formulated  in  Section  II 
for  the  case  of  planar  geometry  and  plane  wave  incidence.  These 
equations,  subject  to  appropriate  boundary  conditions,  are  solved  in 
Section  III.  Therein,  Chebyshev  polynomials  are  used  as  basis 
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functions  for  series  expansions  of  the  diffuse  or  incoherent  inten¬ 
sity.  This  choice  of  basis  functions  reduced  two  of  the  three  linear 
systems  of  equations  which  had  to  be  solved  to  systems  of  two  equa¬ 
tions  each.  Legendre  polynomials  have  been  used  in  a  similar 
fashion  previously  [ZScl.  The  unknown  expansion  coefficients,  how¬ 
ever,  are  related  by  continued  fractions,  whereas  the  expansion  co¬ 
efficient  associated  with  the  Chebyshev  polynomials  are  related  more 
simply  (see  Equations  (25)  and  (40) ) .  In  Section  IV  power  flux  den¬ 
sity  is  defined  in  terms  of  intensity  and  is  shown  to  reduce  to  trun¬ 
cated  series  expansions.  Numerical  results  for  intensity  and  flux 
are  presented  in  Section  V. 

II.  FORMULATION 

A  forest  is  modeled  as  a  slab  region  (Fig.  1)  confined  within 
infinite  planes  z  =  0  and  z  =  d.  Its  vegetation  is  assumed  to  con¬ 
sist  of  a  random  distribution  of  particles  which  scatter  and  absorb 
the  incident  radiation.  The  vegetation  medium  is  assumed  to  be 
statistically  homogeneous  and  to  be  characterized  by  absorption  (cr  ) 

a 

and  scatter  (a  )  cross  sections  per  unit  volume  whose  dimensions  are 

-1  s 

meters  .  Furthermore,  we  assume  that  the  forest  medium  has  an  iso¬ 
tropic  scatter  characteristic,  i.e.,  the  average  scatter  current  will 
radiate  an  equal  amount  of  power  into  each  spatial  direction.  A  uni- 

A 

form  plane  wave  is  assumed  to  enter  the  forest  from  the  direction  ftQ 
defined  by  angles  0Q  and  <j>q  of  a  spherical  coordinate  system  (Fig.  1). 
The  magnitude  of  the  incident  Poynting  vector  or  power  flux  density 
is  given  as  SQ  in  watts  per  meters  squared  and  frequency  is  specified 
by  vQ  in  Hertz. 

As  energy  enters  the  forest  it  is  multiply  scattered  in  all 
directions.  To  describe  this  phenomena  we  use  the  transport  equa¬ 
tion  (also  called  the  equation  of  radiative  transfer) .  The  transport 
equation  is  essentially  a  statement  of  conservation  of  energy  and  has 
as  its  dependent  variable  the  quantity  called  spectral  specific  in¬ 
tensity  1^  (r;ft).  It  is  written  in  this  manner  to  indicate  that  it 
is  both  a  function  of  position  r  =  (x,y,z)  and  of  direction  as  rep- 

A 

resented  by  the  unit  vector  ft;  the  subscript  v  denotes  its  spectral 
nature.  To  define  spectral  specific  intensity  consider  an  element 
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of  area  dA  (Fig.  2)  whose  orientation  is  defined  by  the  unit  normal 
vector  fi.  Assume  that  radiant  power,  i.e.,  energy  per  unit  time  dP 
flows  in  the  direction  specified  by  5  through  the  projected  element 
of  area  dA^  =  dA  cos  6.  The  spectral  specific  intensity  at  the  point 
r  in  the  direction  6  is  then  defined  by  the  relationship  I  (r;ft)  = 
dP/dA^dvdft,  i.e.,  to  be  the  power  per  unit  area  in  watts  per  meters 
squared,  per  frequency  interval  dv  in  Hertz  and  per  solid  angle  dft  in 

A 

steradians  which  flows  in  the  direction  ft  normal  to  the  projected 
element  of  area  dA^. 

Since  the  planar  geometry  and  incident  plane  wave  problem  de¬ 
picted  in  Fig.  1  involves  a  medium  (the  forest)  which  scatters  radi¬ 
ation  isotropically,  the  spectral  specific  intensity  depends  only  on 
the  coordinate  z  and  the  directions  8  and  <|>  .  Since  the  medium  is 
linear  the  scatter  field  has  the  same  frequency  as  the  incident 
radiation.  Hence^ 


Iv(r;  ft)  =  I  (z ;0,<|>)  6  (v-vQ)  , 

and  the  transport  equation  becomes  [OZ,  CZ,  IS1,  Ch] 
j-r  /  2rr  r  IT 


0  al  +  °tI  “  7?  f  Kz;®’,*')  sin  0’  d0'  d<f>',  0<z<d,  (2a) 

'  /n  Jn 


where 

°t  =  °a  +  °s  (2b) 

is  the  total  cross-section  or  extinction  cross-section  per  unit  volume. 
Introducing  the  normalized  quantities 


and  albedo 


T  =  0tZ,  p  =  COS  0 


W  =  as/ot 


allows  us  to  rewrite  Equation  (2)  as 

r 2tt  f  +  1 


where  t 


I  (  t  ;  y  ^ )  dp ' d$ 1 ,  0<t<tq, 


0  -1 


■k)nly  the  coherent  component  of  the  intensity  depends  on  <p ;  the  inco¬ 
herent  component  depends  on  z  and  0;  see  Equations  (13)  and  (14) . 
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For  convenience  we  separate  the  intensity  I  (T;y,<J>)  into  a  sum 
of  two  terms  (see  Fig.  3  and  [Isl]) 


I 


rd  +  rri 


(5) 


After  substituting  Equation  (5)  into  Equation  (4)  1^  is  defined  to 
satisfy  the  relationship 


V 


r2*  I1 
0  -1 


<!d  +  :ri> 


dy ' 


d<j> ' 


(6) 


and  I  ^  to  satisfy  the  simpler  relationship 

=  0. 


dlri  A  T 
y  ~  +  I 


dx 


n 


(7) 


Integrating  Equation  (7)  gives 


(8) 


Comparing  Equation  (7)  with  Equation  (4)  we  see  that  they  would 
be  identical  if  the  integral  terra  were  zero.  The  integral  term 
accounts  for  the  multiscatter  effects  of  the  medium;  without  it  the 
transport  equation  would  simply  describe  the  attenuation  of  energy  as 
it  flows  in  a  particular  direction  through  the  medium  due  to  absorp¬ 
tion  (a  )  and  scattering  (a  )  by  particles  along  a  given  ray  path. 
This  effect  is  described  by  the  solution  given  in  Equation  (8) . 

Hence,  we  call  I  ^  the  reduced  incident  or  coherent  intensity  and  we 
designate  1^,  defined  by  Equation  (6) ,  the  diffuse  or  incoherent  in¬ 
tensity  since  it  accounts  for  multiscatter  effects. 


To  specify  boundary  conditions  required  by  the  slab  geometry  of 
Figure  1  it  is  convenient  to  separate  the  diffuse  intensity  into  two 
components,  a  forward  diffuse  intensity  I^+  and  a  backward  diffuse 
intensity  1^  ,  where 


I 


d 


Jd+  0<6<  7T/2 

.  *  -g  <  0  <  TT 


or  0<y<l 
or  -l<y<0 


(9a) 

(9b) 
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Definition  (9)  is  further  clarified  by  Figure  3  wherein  the  forward 
diffuse  constituent  I^+  is  confined  to  angles  of  scatter  8  between 
zero  and  ninety  degrees  (representing  energy  transport  in  the  positive 
x  direction)  while  the  backward  diffuse  intensity  1^  ranges  over 
ninety  to  one  hundred  and  eighty  degrees  (corresponding  to  energy  flow 
in  the  negative  t  direction) .  Since  these  components  are  generated 
within  the  forest  we  must  have 

I,+  *  0  at  t  —  0  (10a) 

and 

Id"  -  0  at  t  =  tq(  (10b) 

i.e.,  diffuse  intensity  cannot  enter  the  slab  from  the  exterior 
regions . 

Recall  that  we  have  assumed  that  a  plane  wave  enters  the  forest 
in  direction  with  a  power  density  or  Poynting  vector  magnitude  of 
Sq  watts  per  meters  squared.  The  incident  intensity  in  watts  per 
meters  squared  per  radians  squared  is  then 

io  =  sq  6(fi  -  no)  =  so  6(e-eo)S(4>-4>o)  -  so6(v-po)6(.<m.0), 

sin  0  (11) 

where  6(x)  represents  a  Dirac  delta  function  and  ( 0Q ,  4>q )  are  the  di¬ 
rection  angles  of  the  incident  plane  wave.  From  Equations  (5),  (9) 

and  (10)  we  see  that  intensity  in  the  forward  direction  satisfies  the 
continuity  condition 


I  -  I.+  +1=1.  at  i  =  0 
o  d  ri  rx 


(12) 


Using  Equations  (8),  (11)  and  (12)  we  obtain  for  the  reduced  incident 

intensity 

Iri(T;p,*)  =  I0e'T/M=  S0  ■5Cli-P0)‘5(<f’-'t0)e'T/‘,o.  (13) 

Substituting  Equation  n.3)  into  Equation  (6)  and  integrating 
over  <j>  yields  transport  equations  for  the  diffuse  intensity 
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dy ' 


dy  n  + 


^fo 

4tt 


(14) 


0<y<  1. 
-1  0 


The  superscript  designations  (+)  identify  the  appropriate  ranges  of  y 
as  defined  in  Equation  (9) .  The  two  linear  integrodif ferential 
Equations  (14)  subject  to  the  two  boundary  conditions.  Equation  (10), 
constitute  a  complete  mathematical  statement  of  the  problem  which  will 
be  treated  in  the  succeeding  sections. 


III.  REDUCTION  OF  THE  TRANSPORT  EQUATION  TO  SYSTEMS  OF  LINEAR 
EQUATIONS 


The  Transport  equations  (14)  are  solved  by  representing  the 
diffuse  intensities  as  truncated  series  expansions  of  half-range 


Chebyshev  polynomials 


I 


(t;y) 


B^U)  T  “(y  *  »  0  £Vl£l> 
n  n  -10 


(15a) 


with 


Tnt(u)  =  Tn(2y  +  1) 


(15b) 


denoting  two  Chebyshev  polynomials  of  order  n  which  are  complete 
orthogonal  functions  over  the  y  domains  (0,1)  and  (-1,0),  respective¬ 
ly.  Substituting  Equation  (15)  into  Equation  (14)  gives 


N  dB-(i) 

z:s-£ 


<n 

+  l  [B*(t )+B~(t  ) ] [—g— 

n=0 ,2,4,.  .  .  n  _1 


+  *§  -x/y, 


4w 


o  0<y< 1 , 
'-l  0 


(16) 


where  use  was  made  of  the  integral  evaluations 
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[  T*(y)  dy  =  [  T~(y)dy  =  i  f  Tn(x)dx  =  I  n^-l  ’  “  °,2’4,‘" 
0  -1  _1  l  0  ,  11=1,3,5,  .  .. 


(17) 


Using  the  recursive  formulas 


2(2y+l)Tn(2y+l)  =  Tn+1(2y+l)+Tn_1(2y+l)  (18) 


to  replace  the  terms  yTn(2y+l)  in  Equation  (16)  yields 


£  N 

1  itTi+l+Tn-I  J  +  <^2  I?"  +  ^  =2W  n-0,2,4. J  ^ 
n=0 


(19) 


+  5t»  e"T/po  >  0<  y  <1 
-1  0 


TT 


The  orthogonality  properties  of  the  half-range  Chebyshev  poly¬ 
nomials  T*(y)  =  Tn(2y+1)  provide  the  means  to  obtain  the  following 
system  of  2n+2  coupled  first  order  differential  equations  for  the  in¬ 
dependent  expansion  coefficients  Bn(x): 


dB“ 


f-  ±2  T1+4B0  +  ^  ^  cj£  +  C5^T>-  -ir  e'T/li° 

n=o y ^ f 4 1  •  •  • 


WS 


dt 


0,  n=0  ,  (20a) 


dE5  dBr 

2  ~3T~  +  3t  l2 


dB^  + 

-af-  +4BI  -  0 


n  =  1, 


(20b) 


dB 


n-1 


dt 


dB*.  dB; 
n+i  +  2-  n 


di  dx 


+  4Bn  =  0,  n  =  2,3,...,  N-1, 


dB 


N-1 


dx  - 


dBN*  ± 

+  2  +  4B«  =  0,  n=N 


sr 


(20c) 


(20d) 
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In  deducing  the  above  use  was  made  of  the  fact  that  T_^  =  and 

that  the  truncated  series  representations  Equation  (15)  require  the 

coefficients  B  (x)  to  be  zero  for  n>N;  the  latter  was  used  in  writing 
n 

Equation  (20d) . 

Orthogonality  of  the  Tn(2y+1)  functions  is  also  used  to  obtain 
from  Equations  (10)  and  (15)  the  boundary  conditions 


Bn(0)  =  0 

vv  -  °- 


for  all  n 


(21a) 

(21b) 


The  total  solution  to  Equation  (20)  subject  to  boundary  con¬ 
ditions  Equation  (21)  consists  of  homogeneous  and  particular  solutions 

B~  (x)  and  B-^3  (x),  respectively,  i.e., 
n  n 


-  n±b 


Bn(T)  =  Bn  (t)  +  Bn  (t)’  n  = 


(22) 


The  coefficients  B~  (x)  satisfy  the  linear  system  of  Equations  (20) 
with  Sq  set  equal  to  zero  while  (x)  satisfy  the  inhomogeneous  sys¬ 

tem  (SQ^0)  and  have  the  same  functional  dependence  on  x  as  the  forc¬ 
ing  function  (see  Equation  (42b)). 

A.  Homogeneous  Solution  to  Equation  (20)  (SQ  =  0) 

Since  the  system  of  equations (20)  involves  only  first  order 
derivatives  j_n  T  f  homogeneous  solutions  must  depend  exponentially  on 
x.  Hence,  it  may  be  assumed  that,  in  general,  the  homogeneous  solu¬ 
tion  to  Equation  (20)  takes  the  form 


B‘h(T)  .[,*  e~T/sk  , 
n  n,k  ' 

k=0 


(23) 


where  a~  ,  and  s,  are  unknown  constants  and  integer  k  accounts  for  all 
n,k  k 

possible  solutions.  Substituting  Equation  (23)  into  Equations  (20) 
and  (21)  gives  the  linear  system  of  equations 
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<N 


alk-  2(2sk+1)a0,k~2sk  W  ^o,2,. 


2a0,k2(2Sk"i)al,k  ^k  °’  n=1  * 


aiLk=  o,  n=0  ,  (24a) 
n-l 


(24b) 


an-l,k  "2(2sk+1)an,k  +  an+l,k 


=  o,  n=2 , ...»  N-l »  (24c) 


and 


4-1, k  -2(2sk?l)a*  k=°,  n-H 


(  24d ) 


Using  the  definitions  of  the  Chebyshev  polynomials  of  the 
second  kind  [ASt] ,  the  three  Equations  (24b) ,  (24c)  and  (24d)  yield 

expressions  for  all  coefficients  in  terms  of  the  two  coefficients 

aN,k: 


an,k  =  DN-n<2V1)aN,k  ■  ■>=<>. I.2 . N. 

1,  n=0 

2,  n+0  . 


en  =< 


(25a) 


(25b) 


Substituting  Equation  (25)  into  Equation  (24a)  and  utilizing  the  re- 
lationships 


Un+l(x)  =  2x  U-(x)  "  U— i(x> 


n 


n-l 


T„.,(x)  =  x  UN(x)  -  UM_,(x) 


N+l 


N-l 


(26a) 

(26b) 


result  in  the  simpler  system  of  equations 
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A(sk)ak  =  0, 


where 


~ (sk>  “ 


al(sk)  ao(0 


2k' 


;(sk) 


-k  =  \  a 


(27a) 


(27b) 


sk  =  2sfc  +1,  a2(x)  =  a(x)-^WB(x),  a2(x)=skW6(x),  (27c) 


£Nu  (x) 

a(x)  =  Tn+1(x),  8(x)  =  -  l^(x)  +  2S^(x),Se(x)=£j-±— ;pCi  (27d) 

n-2 , 4 , • « • 

For  non-trivial  solutions,  det  A(Sk)  =  0;  yielding  the  charac¬ 


teristic  equation 


TN+i<2V1)W2sk-1>+V'  tWV” fw^o'vl 

f  -1  (28) 

+  TN+1(2sk-l)|-U(2s+l)+2Se(2sk+l)|]  =0- 

2 

This  eigenvalue  equation  is  a  polynomial  in  sk  with  real  coefficients 

of  order  N+l  and  provides  2N+2  non-zero,  distinct,  real  eigenvalues 

s,  =  +  s  ,  s  >0,  m=0,.  .  .  ,N .  Physically,  real  eigenvalues  are  to 
Jc  —  m  in 

be  expected  since  the  associated  eigenfunctions  are  orthogonal  and 
their  sum  represents  specific  intensity  or  real  power.  Combining 
Equation  (28)  with  one  of  the  equations  of  Equation  (27)  yields  for 
each  eigenvalue  sk : 


iSJt  .  TN+l(2sk-I)  ,  k  -  ±m,  m=0 . N. 

Vk  W2Sk+1) 

For  k  =  m,  i.e.,  for  positive  eigenvalues  sk  =  sm>0  assume  that 

a*  =  am  Txjj_i  (2s  +  1) ,  a  =  constant. 

N,m  m  N+l  m  m 


(29) 


(30) 


It  then  follows  using  Equations  (29)  and  (30)  that 


aw  ™  =  a~  T  ( 2s  -1). 
N,m  m  N+lv  m  7 


(31) 
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(32) 


Combining  Equations  (30)  with  (31)  gives  lor  s.=s  >0 

j\  in 


a„  m  =  a  T._l.(2s  ±1). 
N,m  m  N+l  m  7 


Similarly,  for  k  =  -m,  i.e.,  for  negative  eigenvalues 

s,  =  s  -  -  s  it  can  be  shown  that 
k  -m  m 

a*  =  a  m  T  (2s<+l),  a  m  *  constant. 
N,— m  —  m  N+l  m  — ro 


(33) 


Using  Equations  (32)  and  (33)  with  (25)  permits  Equation  (23)  to  be 
written  in  terms  of  positive  eigenvalues  only: 

Bf<T>  *  5  en  £  tam  V'V)  "h<2V')e'’,S" 


+<-l  >%  W2sm?1) 


In  Equation  (34)  we  have  used  the  identity  U  (-x)  =  (-1)  U  (x)  and 

n 

have  set  b^  =  (-1)  a_m«  Hence,  with  Equation  (15)  the  most  general 
homogeneous  solution  to  Equation  (14)  is: 


lT  Bnh(T)  Tn(2^1) 


(35) 


with  (t)  given  by  Equation  (34) . 

i  y. 

Simpler  expressions  for  are  obtained  by  substituting 

Equation  (34)  into  (35) ,  interchanging  the  order  of  the  summation  op¬ 
erations  and  using  a  closed  form  expression  to  replace  the  series 
which  is  summed  over  the  integer  n.  The  result  is 


I  t  ame-T/s">  T  (2s  ±1)  [  T^l(2S",'l'1)~TW-I(2>''fl) 

4  m-0  W2*.*”!  Vp  ] 

_  T  (2s  ±1)+(-1)N  T(J+1(2yIl) 

•c  b  e  T/S  T  (2s  +1)  C  iSti - 5! - ESI -  ] 

m  '  m  an+1  m  '  sm-u  J 

(36) 
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B.  Particular  Solution 


To  obtain  the  particular  solution  we  assume  that  B^(t)  has  the 
same  x-  dependence  as  the  forcing  term  in  Equation  (20a) ,  namely. 


B;p(t)  =  A*  e-T/V 
n  n 


(37) 


Substituting  into  Equation  (20)  yields  the  following  linear  system 
of  equations  for  the  coefficients  An: 


N  a++  A- 

-2(2  po+l)  A±  -2^W2  +  Po^=  0,  n  =  0, 

n  -1 

n=0,2,4  ,  •  •  • 


(38a) 


2A^  -  2(2po+])A“  +  A^  -  0 


n  =  1  , 


(38b) 


A^  -2C2y0-fl)A^  +  A*+1  =  0  ,  n  =  2,  . . . ,  N-l  ,  (38c) 


and 


Aj^  -2(2yQTl)Aj  =  0  ,  n  =  N 


(38d) 


Since  Equations  (38)  are  identical  in  form  to  Equations  (24)  except 
for  the  forcing  term  yoSoW/7r  they  can  be  solved  in  like  fashion  to 
yield  the  system  of  two  equations  for  the  two  unknowns  aN: 

i-C*)  •  1-  <  l  ).  (39) 

where 

Ad  "  i  en  °K-n(2't>:fl,AJ  .  n  -  0,1,2, ...»  ,  (40) 

/it 

with  defined  in  Equation  (25b) . 

The  matrix  A  (P  )  is  identical  to  A(s,  )  in  Equation  (27b)  but  with 
s,  replaced  by  M  and  s;~  by  pi  —  =  2 pi  +1.  The  solution  to  Equation 
(39)  involves  a  simple  2x2  matrix  inversion  from  which  it  is  found 
that 
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(41) 


AN  - 


ws. 


TT 


W2%^> 

det  A(uq) 


where  det  A(yQ)  is  identical  to  the  left  hand  side  of  Equation  (28) 
but  with  s,  replaced  by  y  . 

K  O 


Combining  Equations  (37)  and  (39)  with  Equation  (41)  gives  the 
particular  solution 


with 


B^(t)  Tn(2y+1) 

n=0 


(42a) 


(T) 


y  WS 
.  Ho _ o_ 

'n  2tt 


TN+l(2yoi;L) 


det 


lW2V1>e"T/’J° 


(42b) 


As  was  done  for  the  homogeneous  solution,  the  summation  over  n  in 
Equation  (42a)  can  be  replaced  by  a  closed  form  expression  which  re¬ 
sults  in 

,±P_  4  e~T/v->  w2'Vl>  -  w211*”  <43) 

Id  "  4  I-  yQ  ~y  J 

with  A^-  given  by  Equation  (41)  . 

C .  Total  Solution 

The  complete  solution  for  the  diffuse  intensity  is  the  super¬ 
position  of  the  homogeneous  and  particular  solutions,  given  either  by 
Equations  (35)  and  (42)  or  (36)  and  (43) .  The  unknown  coefficients 
am  and  bm  are  determined  by  application  of  the  boundary  conditions  in 
Equation  (21)  which  we  write  as 


B**1  +  B+P  =  0  at  t  =  0  for  0  <  y  <  1 
n  n  -  - 


(44a) 


and 


=  t  for  -1  <  y  <  0  . 


0  at  t 


(44b) 


Substitution  of  Equations  (34)  and  (42b)  into  Equation  (44)  and 
equating  coefficients  of  Tn  (2y+l)  to  zero  (these  polynomials  are 
complete  for  the  y-ranges  involved)  provide  the  following  system  of 
inhomogeneous  linear  equations  from  which  to  determine  the  2N+2  un¬ 
known  constants  a  and  b  : 

m  m 

EfVW 'V1 >°N-n ( 2sn,  >Vl ' <  ^'W2*^ 

m=0  , 

<45a> 

E  IVn+1<2s  ^UN.n(2sm-i)eVsn, 

rn=0 

=  -A"n  e  To/Uo  UN_n(2yQ+l)  (45b) 

with  given  by  Equation  (41)  and  n  =  0,...,N. 

Although  we  have  formally  addressed  only  the  slab  problem,  the 
above  analysis  applies  in  a  simple  fashion  to  *-he  half  space  problem 
(To  =  °°  )  .  From  Equation  (45b)  we  see  that  for  all  coefficients 

b^  must  be  identically  zero  in  order  that  the  system  of  equations  re¬ 
main  bounded.  Equation  (45a)  alone  yields  enough  relationships  from 

which  to  determine  the  remaining  nonzero  coefficients  a  . 

J  m 

IV.  FLUX  DENSITY 

A  second  fundamental  quantity  used  in  transport  theory  to  de¬ 
scribe  power  flow  properties  is  the  flux  density.  With  reference  to 
Figure  2  and  the  definition  of  specific  intensity,  the  forward  flux 
density  F  +  (r;fi)  is  defined  [Isl]  as  the  amount  of  power  that  flows 
through  the  area  element  dA  in  all  directions  0  <  8  <  tt/2  where  0  is 
measured  from  the  unit  vector  n  normal  to  the  area  element,  i.e., 

C  2tt  tt/2 

F+(r;n)  =  I(r;fi)ft.n  dQ  ,  fl.n  =  cos  6.  (46a) 

■  0  -*0 

Flux  is  in  units  of  watts  per  meters  squared  and  solid  angle  dft  = 
sin  0d0d<}>.  In  like  fashion,  we  also  define  a  backward  flux  density 
F~(r;fi)  in  terms  of  power  which  flows  through  dA  but  in  the  -n  di¬ 
rection,  i.e.  , 
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2TT  f  TT 

F”(r;n)  =  I(-r;6)n.  (-n)dfi  . 

■  0  ■*  ff/2 


(46b) 


For  the  slab  geometry  and  plane  wave  incidence  of  Figure  1,  the 
above  flux  densities  depend  only  on  z  (or  t  )  and  can  be  written 
using  Equations  (8)  and  (15)  as 


F+(t)  =  P^CxJ+F+tT),  r'(T)  = 


(47a) 


where 


F^t) 


2u  ,1 


yl  ,  (T;<j>,u)dyd<t>=Se  T^o  6(<p-<p  )d<J>  6(vi-p  )ydy 
x  J*  o  o  o 


o  o 


=  y  Se 
o  o 


— x /Pc 


(47b) 


F^(t)  * 


■2tt  cl  N  ,1 

pI*(T;y)dyd<fr  =£  B*(t)  [  2tt  yT^  (2y-l)dy], 

n  A  n=0 

0  0  n 


(47c) 


F-(t) 


■2tt  cO  N  0 

Vld(r,v)dud4  =  -]T  Bn(t)[  271  yTn(2y+l)dyJ , 

•  n=0  J_5 

0-1 


(47d) 


where  y  has  been  substituted  for  9 . 

Evaluations  of  the  above  integrals  can  be  shown  to  yield 


PTn 


(2y-l)dy  * 


n-1 

-  TT 

1^4 


,  n  =  0,2,4, 


,  n  »  1,3,5, . 


yT  (2y  +1  )dy 
-1  n 


,  n  =  0,2,4, ..  . 
,  n  —  1,3,5,... 
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Hence 


F±(t) 

where  the  summations  terminate  at  N  or  N-l ,  N  being  either  even  or 
odd  and  B^(x)  is  given  by  Equation  (22)  with  available  from 

Equation  (34)  and  (  x  )  from  Equation  (42b)  . 

A  convenient  way  to  express  the  above  result  after  interchanging 
the  order  of  the  summation  is  to  use  the  quantities  defined  in 
Equation  (27d) ,  namely  a(x),  3  (x)  and  Se  (x)  along  with  the  odd  sum¬ 
mation  term 

<N  9 

SQ(X)  -£  (n  -4)  UN_n(x)  (51) 

n=l }3}5>a  *  • 


=  -  7T 


<N  <N 

£  “K  Ei(J)  'jr4  ^(T) 

n=0,2,4,n.  ..  n=l,3,5,  ... 


(50) 


so  that  Equation  (50)  becomes 


N 


F~(T>  ■  ♦  2  So(x^))e-T/sm 


+  bma(xj){8(x*)  +2  SD(x^)}e  T/sm  ] 


m ' 


(52a) 


U  WS 


•  2°det  A ~  “(VoHSCy*)  +  2  So(y*)}e-T/',o, 


where 


2s  +1, 
m  ’ 


+1- 


(52b) 


V.  NUMERICAL  RESULTS 

Based  on  the  Chebyshev  method,  data  is  obtained  for  the  diffuse 
intensities  and  the  power  fluxes  for  both  the  half-space  (x  =  °°)  and 
slab  configurations.  Selected  numerical  results  for  flux  incorporate 
a  comparison  with  calculations  using  first  order  theory  expressions 
(see  ( Isl ] )  . 


As  can  be  seen  from  Equation  (4)  cr  {  1  •: .  the  forest  is  charac¬ 
terized  by  the  single  parameter  W  called  the  albedo,  which  is  defined 
in  Equation  (3b).  In  general,  it  can  vary  from  zero  (no  scattering) 
to  unity  (no  absorption) .  We  have  chosen  W  to  range  from  a  low  of 
0.25  to  a  high  of  0.95;  the  former  describes  a  medium  which  mostly 
absorbs  energy  whereas  the  latter  characterizes  a  predominantly 
scattering  environment.  A  plane  wave  is  assumed  to  enter  the  forest 
at  incidence  angles  of  6q  =  0°,  30°  and  60°  for  the  half-space  geo¬ 
metry  and  at  0  =0°  for  the  slab  geometry.  The  normalized  slab 

thickness  xq  varies  from  0.2  to  20.  All  intensity  curves  are  normal¬ 
ized  by  setting  Sq/4it  equal  to  unity  while  all  power  flux  curves  are 

normalized  with  S  set  to  unity.  These  normalizations  were  chosen  to 

o 

be  compatible  with  the  units  of  specific  intensity  which  is  power  per 
unit  area,  per  unit  solid  angle  and  of  power  flux  which  is  power  per 
unit  area. 

In  Figures  4  and  5,  the  forward  I*  (x;  0°)  and  backward  1^ 

(t;180°)  diffuse  intensities  in  the  0=0°  and  180°  directions,  re¬ 
spectively,  are  plotted  for  the  three  angles  of  incidence  0q  =  0°, 

30°  and  60°  as  a  function  of  penetration  depth  x  for  albedos  W  =  0.25, 
0.50,  0.75,  0.90  and  0.95.  For  very  small  x  values  the  effect  of 
multiscattering  is  small.  The  flow  of  energy  is  mainly  confined  to 
the  coherent  (reduced  incident)  intensity  (not  shown) .  Soon,  however, 
multiscattering  of  the  energy  takes  place.  Energy  transfers  into 
the  diffuse  or  incoherent  state.  It  quickly  builds  up  to  a  maximum, 
having  started  from  zero,  but  then  decreases  with  x  due  to  absorption 
loss  along  the  longer  multiscatter  path  lengths  traversed  in  going 
from  x  =  0  to  an  observation  point  x.  Meanwhile,  the  coherent  energy 
decreases  more  rapidly  (see  e.g.  curve  5  in  Fig.  15. c)  since  it  loses 
energy  both  due  to  scattering  out  of  the  direction  of  flow  as  well  as 
to  absorption  losses  along  the  prescribed  path.  Furthermore ,  the 
maximum  forward  diffuse  intensity  is  larger  and  occurs  deeper  inside 
the  half-space  region  for  stronger  scattering  (larger  W)  media.  It 
also  takes  longer  for  intensity  to  attenuate  and  its  rate  of  atten¬ 
uation  is  smaller  when  W  is  larger. 

Observe  further  that  the  forward  diffuse  intensity  in  the 
0=0°  direction  i^  larger  for  smaller  incidence  angles.  This  can 
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be  explained  by  noting  that  at  a  given  value  of  x  more  energy  is 
available  from  the  coherent  intensity  I  ^  to  be  mul tiscattered  when 
0q  is  smaller.  For  larger  x  values,  the  rate  of  decrease  of  1^  (x;0°) 
depends  only  on  the  albedo.  This  follows  from  the  asymptotic  solu¬ 
tion  (see  [ C Z ])  or  from  the  curves  in  Figure  4  for  W  =  0.9  and  0.95 
which  have  constant  slopes  near  x  =  10  (the  asymptotic  x  range  is  not 
reached  for  curves  with  smaller  values  of  W) .  However,  for  smaller  x 
values,  incidence  angle  does  influence  the  rate  of  attenuation.  The 
curves  in  Figure  4  indicate  that  smaller  rates  of  attenuation  are 
associated  with  smaller  incidence  angles.  This  also  follows  from 
Equation  (43)  where  the  exponential  attenuation  (exp ( -x/y  ))  is  a  dom¬ 
inant  decay  factor?  for  y  equal  to  unity,  i.e.,  for  normal  incidence 
(0  =  0°) ,  the  attenuation  rate  in  the  x  direction  is  smaller  than 

when  y  is  smaller, 
o 

Figures  5  depict  the  x  dependence  of  the  normalized  (So/4ir  =  1) 
backward  diffuse  intensity  (in  decibels)  which  flows  in  the  9  =  180° 
direction.  Again  we  see  that  the  larger  the  albedo  the  stronger  the 
incoherent  intensity  which  is  characteristic  of  increased  scattering 
taking  place  in  the  forest.  Note  that  the  strongest  backward  diffuse 
intensity  occurs  at  the  air-forest  interface  x  =  0.  This  was  antic¬ 
ipated  because  only  the  region  to  the  right  of  a  particular  value  of 
x  contributes  to  the  backward  intensity  1^  at  the  value  x?  hence, 
since  every  particle  in  the  entire  half-space  (T>0)  scatters  energy, 
a  portion  of  which  can  be  assumed  to  reach  the  interface,  maximum 
backscatter  intensity  occurs  there.  As  in  the  case  of  forward  scat¬ 
ter,  the  angle  of  incidence  also  affects  the  amount  of  energy 
which  backscatters .  We  see  that  for  larger  6q  less  energy  travels  in 
the  6  =  180°  direction;  this  is  due  to  the  fact  that  a  larger  0Q  re¬ 
sults  in  a  smaller  coherent  intensity  reaching  the  depth  x  in  the 
half-space  region;  hence,  a  smaller  coherent  intensity  penetrates  the 
region  to  the  right  of  x.  With  less  energy  available  from  the  coher¬ 
ent  energy  state  and  thus  for  larger  0q  less  scatter  energy  flows  in 
the  backward  direction. 

Figures  6  through  8  are  plots  of  the  normalized  diffuse  inten¬ 
sity  (in  decibels)  versus  scatter  angle  0  for  albedos  W  =  0.25,  0.50, 
0.75,  0.90  and  0.95  and  for  different  penetration  depths  x.  in 
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Figure  6,  I^  (x; 0)/I+ (x; 90u)  is  plotted  for  the  smaller  depths  (x=0.05, 
0.1,  0.2/  0.5  and  1.0)  while  in  Figure  7,  (x ; 0) /I+ (x ; 0°)  versus 

0  is  plotted  for  the  larger  depths  (x  =  1,  2,  5,  10,  20).  In  Figure 
8,  I(j  (x  >  0) /I(j  (x;  90°)  versus  0  is  plotted  for  the  three  representative 
values  x  =  0.1,  1.0  and  10.  Also  drawn  are  curves  for  different  in¬ 
cidence  angles  (0  =  0°,  30°  and  60°) . 

Consider  first  Figure  6.  It  is  immediately  seen  that  for  small 
penetration  depths,  maximum  intensities  flow  in  directions  near  0  = 
90°.  As  x  increases,  the  maxima  shift  toward  the  0=0°  direction. 
Furthermore,  for  larger  albedos  the  relative  peak  intensity  values 
diminish  and  occur  at  larger  0  values  (compare,  for  example,  curve  5 
in  Figure  6. a  to  curve  5  in  Figure  6.m).  Note  also  in  each  figure, 
and  for  the  larger  x  values  particularly,  that  the  difference  in  in¬ 
tensity,  as  measured  from  one  curve  to  the  next  at  a  given  value  of 
0,  decreases  for  larger  albedos,  i.e.,  the  curves  become  more  com¬ 
pressed  (compare  the  curves  of  Figure  6. a  to  the  curves  of  Figure  6.m 
to  see  this  effect) .  Furthermore,  at  a  given  observation  point  x  we 
see  that  as  0Q  increases  the  relative  peak  intensity  rises  and  shifts 
further  away  from  the  vertical;  for  illustration  refer  to  the  x=l, 
0o=O°  and  60°  curves  of  Figures  6.m  and  6.o,  respectively.  Likewise, 
for  a  given  incidence  angle  0Q,  as  x  increases  this  same  shift  occurs 
as  is  vividly  portrayed  in  the  curves  of  Figures  7. a  through  7.o 
where  the  relative  maximum  are  seen  to  shift  to  the  0=0°  direction. 

Changing  0Q  seems  to  have  a  minor  effect  on  this  shift  in  maxi¬ 
mum  intensity,  but  an  increase  in  x  plays  a  significant  role.  Recall 
that  a  major  portion  of  the  energy  which  scatters  in  a  particular  di¬ 
rection  0  through  the  observation  point  x  comes  from  the  more  direct 
or  shorter  paths  traversed  between  a  scatterer (hit by  the  incident 
radiation)  and  the  observation  point.  Evidently,  for  small  x  and  no 
matter  whether  0q  be  small  or  large,  maximum  intensity  flows  in  di¬ 
rections  near  the  vertical  (0  =  90°)  .  This  results  because  contri¬ 
butions  to  the  energy  flow  through  the  point  x  near  the  air-forest 
interface  in  directions  0Q  are  overwhelmed  by  the  numerous  energy 
contributions  coming  from  the  many  scatterers  in  the  infinite  domain 
of  the  thin  slab  region  between  the  interface  (x  =  0)  and  the  obser¬ 
vation  point  x;  the  bulk  of  the  energy  must  travel  in  directions  near 
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9  =  90°  in  order  to  pass  through  t .  For  fixed  t  and  for  increasing  6q 
more  energy  transfers  into  the  diffuse  intensity  sooner  (i.e.,  at 
smaller  depths  x).  Hence,  more  energy  is  now  available  to  be  scatter¬ 
ed  into  the  0  =>  0°  direction  to  counteract  the  energy  contribution 
which  enhances  vertical  flow.  However,  the  contributions  in  the  nor¬ 
mal  direction  for  small  t  values  are  small  and  the  effect  of  shifting 
the  intensity  away  from  the  vertical  direction  is  not  strongly  felt 
until  x  becomes  larger. 

For  larger  values  of  t,  Figures  7. a  through  7.o  show  that  the 
angular  distribution  of  intensity  evolves  into  a  narrow  beam  of  radi¬ 
ation  directed  in  the  0=0°  direction.  This  can  be  understood  by 
again  considering  a  point  on  the  x-axis,  but  this  time  far  from  the 
air-forest  interface.  The  strongest  scattered  radiation  to  reach  the 
observation  point  x  is  now  seen  to  come  from  the  region  adjacent  to 
the  x-axis.  This  occurs  because  energy  originating  there  suffers  the 
least  loss  (due  to  absorption)  as  paths  leading  to  the  observation 
point  from  this  region  are  smaller  than  from  elsewhere.  In  Figures  7 
it  is  seen  that  for  the  smaller  albedo  cases,  indicative  of  larger 
absorption  losses,  beam  narrowing  occurs  at  shorter  penetration 
depths  when  compared  to  the  larger  albedo  cases.  Hence,  the  forward 
diffuse  intensity  exhibits  focusing  into  the  0=0°  direction  at 
large  x  values  regardless  of  the  direction  of  the  incident  radiation. 

In  Figures  8,  the  relative  backward  diffuse  intensity  I^(9;x)/ 

Id (90°;x)  versus  scatter  angle  0  for  the  half -space  problem  is  seen 
to  display  a  maximum  intensity  flowing  in  the  0  =  90°  direction  (ex¬ 
cept  for  the  x  =  0.1,  0q  =  0°  curves  of  Figures  8.j  (W  =  0.9)  and 
8.m  (W  =  0.95)  and  the  x  =  0.1,  0  =  30°  (W  =  0.95)  curve  of  Figure 
8.n)  and  to  decrease  monotonically  to  a  minimum  value  in  the  0  =  180° 
direction.  Also,  for  the  small  albedo  case  (W  =  0.25)  and  for  in¬ 
cidence  angles  0q  =  0°  and  30°  (Figs.  8. a  and  8.b,  respectively)  the 
relative  backward  diffuse  intensity  is  shown  to  be  essentially  inde¬ 
pendent  of  x;  this  is  also  true  for  the  x  =  0.10  and  1.0,  0  =  60° 

o 

curves  of  Figure  8.c  and  for  the  0q  =  60°  curves  of  Figures  8.f,  8.i, 
8.1  and  8.o. 

The  curves  show  that  as  energy  progresses  through  the  forest  it, 
at  first,  scatters  predominantly  in  the  vertical  direction  and 


v. 
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eventually  proceeds  in  the  forward  normal  aiiection,  losing  energy  all 
the  time.  Backscatter  at  x  comes  from  the  entire  region  to  the  right 
of  x .  In  a  small  layer  to  the  right  of  x,  energy  is  the  strongest  and 
provides  the  most  significant  backscatter  contribution.  This  contri¬ 
bution  flows  predominantly  in  the  near  vertical  direction.  As  0  in¬ 
creases  from  90°  to  180°,  the  backscatter  diffuse  intensity  falls  be¬ 
cause  energy  contributions  in  these  directions  (from  successive  layers 
to  the  right  of  the  initial  layer)  are  relatively  weaker.  Hence,  in 
these  directions,  the  intensity  is  less  than  in  the  vertical  direction 
as  most  figures  show  and  the  x -dependence  is  essentially  that  of  the 
vertical  direction. 

For  small  0  observe  that  the  relative  backward  diffuse  intensity 
o 

at  small  penetration  depths  is  larger  than  at  large  penetration  depths. 
However,  at  larger  incidence  angles  0  the  reverse  occurs;  less  rel¬ 
ative  energy  flows  in  the  0  =  180°  direction  at  smaller  x  values  than 
at  larger  x  values.  At  larger  angles  0q  more  energy  is  lost  from  the 
coherent  intensity  due  to  absorption  and  scatter  over  a  small  distance 
x  than  occurs  at  the  smaller  angles  0  (note  that  if  £  is  the  path 
length,  x  =  £  cos  0q;  hence,  if  x  is  fixed,  £  is  larger  when  0Q  is 
larger) .  Therefore,  less  energy  is  available  to  be  scattered  into 
the  backward  direction  when  0Q  is  large.  For  small  x,  most  of  the 
energy  remains  still  with  the  coherent  intensity  and  the  diffuse  in¬ 
tensity  strongly  scatters  near  the  vertical  direction  regardless  of 
0Q.  Therefore,  for  small  x  we  expect  the  same  for  the  back  diffuse 
intensity.  However,  for  small  angles  considerable  energy  enters  the 
region  to  the  right  of  x;  the  smaller  t,  the  more  energy  enters  the 
larger  region  to  the  right  of  x  and,  therefore,  contributes  more 
scatter  energy  backward.  For  sufficiently  large  x,  the  forward  energy 
flow  becomes  essentially  forward  diffuse  energy.  It  is  a  beam  of  ra¬ 
diation  centered  around  the  x-axis.  The  strongest  backscatter  inten¬ 
sity  now  travels  near  the  x-axis  but  in  the  0  =  180°  as  path  lengths 
are  smallest  in  this  direction.  For  larger  0Q  this  effect  is  reached 
sooner  in  the  forest  than  when  0Q  is  smaller.  Hence,  we  see  a  strong¬ 
er  relative  diffuse  intensity  at  larger  x  values  when  0Q  is  large. 

In  Figures  9  through  12  intensity  curves  are  drawn  for  the  slab 
problem.  All  curves  are  very  similar  to  those  we  have  already 
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discussed  for  the  half-space  problem.  Note  in  Figure  10  that  1^  goes 
to  zero  as  t  approaches  xq  as  boundary  condition.  Equation  (10b),  re¬ 
quires. 

Power  flux  density  is  plotted  versus  x  in  Figures  13  through  20 
(using  the  Chebyshev  method) .  Curves  for  the  half-space  problem  are 
given  in  Figure  13;  the  remaining  figures  depict  results  for  the  slab 
problem.  Values  for  albedo  range  from  0.5  to  0.95.  A  normalization 
of  incident  power  flux  Sq  equal  to  unity  is  chosen  for  all  flux 
curves.  In  Figures  13-17,  normal  incidence  (0  =  0°)  is  chosen.  In 

the  remaining  three  flux  figures.  Figures  18-20,  0q  is  taken  to  be 
forty-five  degrees. 

In  Figure  13. a,  the  normalized  forward  diffuse  flux  is  seen  to 
quickly  reach  a  maximum  and  then  to  decay  exponentially  with  x.  The 
flux  is  larger  and  its  rate  of  decay  with  x  slower  when  the  albedo  is 
bigger.  This  is  also  true  for  the  backward  diffuse  density  curves  in 
Figure  13. b.  Such  behavior  is  characteristic  of  flux  density  curves 
for  the  slab  problem  as  well,  but  a  dependence  on  xq  enters  the  prob¬ 
lem.  In  particular,  for  small  x  the  flux  F+  monotonically  increases 
from  zero  to  a  peak  value  at  xq  (see  Figs.  14. a  and  15. a  for  xq  = 

0.25  and  1.0,  respectively).  For  larger  values  of  xq  (see  Fig.  16. a 
(xQ  =  5)  and  Figure  17. a  (x  =  20)),  a  maximum  is  reached  and  expo¬ 
nential  decay  then  follows.  For  the  slab  case,  the  backward  diffuse 
flux  density  F^  monotonically  increases  from  zero  at  xq  to  a  peak 
value  at  the  air-forest  interface  (x  =  0) .  This  latter  effect  is 
reasonable  since  it  shows  the  flux  density  to  be  a  maximum  in  the 
backward  direction  when  the  entire  forest  region  is  involved. 

For  the  normal  incidence  case  we  see  in  Figures  14. c,  15. c,  16. c, 
and  17. c  the  exponential  attenuation  of  both  the  coherent  flux  den¬ 
sities  F^  and  total  forward  flux  densitities  Fr^  +  F^  for  different 
albedos.  The  former  (Fr^)  curve  5  in  all  the  figures  exponentially 
attenuates  the  most  rapidly.  This  occurs  because  the  coherent  flux 
density  always  loses  energy  due  to  absorption  and  out-scattering 
whereas  the  incoherent  flux  density  F+ ,  which  also  loses  energy  due 
to  these  effects,  gains  energy  due  to  in-scattering  and  hence  does 
not  attenuate  as  rapidly  as  F^.  Note  that  the  results  for  the  thick 
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slab  tq  =  20)  drawn  in  Figure  1/  are  very  similar  to  the  results  for 
the  half-space  problem  presented  in  Figure  13;  of  course,  the  back¬ 
ward  diffuse  flux  density  F^  is  zero  at  x  =  xQ  for  the  slab  problem 
whereas  it  monotonically  goes  to  zero  as  t  ■+  °°  for  the  half-space 
problem. 

For  forty-five  degree  incidence  angles,  it  is  seen  from  Figure 
18  that  for  small  t  the  coherent  flux  density  is  almost  constant 
over  the  slab  width.  In  addition,  forward  and  backward  flux  densities 
are  nearly  identical  for  the  two  different  incidence  angle  cases  0q  = 

0°  and  45°  (compare  curves  of  Figure  14  to  those  of  Figure  18) . 

Hence,  for  thin  slabs  the  incidence  angle  does  not  affect  flux  den¬ 
sity  very  much.  For  larger  vales  of  tq  (compare  Figure  17  to  Figure 
20)  the  effect  of  changing  6q  is  noticeable ,  but  again,  not  signifi¬ 
cantly.  Observe  further  that  this  effect  is  more  pronounced  when  W 
is  smaller.  This  is  also  true  for  intensity  as  can  be  seen  by  com¬ 
paring  Figure  7. a  with  7.e.  Recall  that  a  smaller  W  means  less  scat¬ 
tering  and  more  abso’-  t  .on.  F*  and  Fd  are  then  smaller  at  larger  0Q 

since  path  length  i-  larger  and  attenuation  stronger  when  0q  is 
larger. 

Flux  curves  are  drawn  in  Figure  21  from  calculations  based  on 
the  first  order  theory  [Isl] .  In  this  theory  the  forward  flux  is  in¬ 
dependent  of  slab  width  i  .  Hence,  Figures  21. a  and  21. b  remain  valid 
for  both  the  half-space  and  slab  problems.  However,  slab  thickness 
does  affect  the  calculation  of  backward  flux  density.  This  is  evi¬ 
dent  by  comparing  Figure  21.c  to  21.d  and  21.e.  The  former  shows 
F^  curves  tor  the  half-space  problem  with  W  =  0.5,  0.75,  0.90,  0.95 
and  0q  =  0°  whereas  the  latter  depicts  backward  flux  curves  for  two 
slabs  (xQ  =  1  and  20,  respectively)  with  W  =  0.5,  0.9  and  0Q  =  0°. 
Figure  21. f  depicts  the  coherent  flux  density  (curve  5)  as  well  as 
the  total  forward  flux  densities  as  a  function  of  t,  valid  for  both 
the  half-space  and  slab  configurations.  Evidently,  F  ^  dominates. 

The  general  characteristics  of  these  first  order  flux  curves  agree 
with  the  more  accurate  ones  presented  previously  using  the  Chebyshev 
method  and  give  additional  support  to  the  latter  theory.  Note  that 
the  Chebyshev  method  yields  a  more  distinctive  effect  of  varying  al¬ 
bedo  values  than  does  first  order  theory  (as  can  be  seen  by  comparing 
Figure  17  to  Figure  21) . 
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In  Tables  1,  2  and  3,  all  positive  eigenvalues  of  the  character¬ 
istic  equation  are  listed.  They  occur  between  0  and  1  except  for  the 
largest  one  which  is  larger  than  unity.  For  larger  W,  the  largest 
eigenvalue  is  closer  to  the  root  of  the  asymptotic  eigenvalue 
equation,  i.e.,  tanh  (1/Ws)  =  1/s  (see  [C Z]  or  [ Oz ]).  The  larger  the 
value  of  W,  the  larger  the  value  of  the  largest  eigenvalue.  The 
largest  eigenvalue  approaches  infinity  when  W  is  equal  to  unity. 

The  smallest  eigenvalue  approaches  zero  as  N  approaches  infinity. 

VI.  CONCLUSIONS 

Diffuse  intensity  curves  in  the  forward  and  backward  directions 
were  presented  for  various  values  of  albedo  (from  a  low  of  0.25  to  a 
high  of  0.95).  A  plane  wave  was  assumed  to  enter  the  forest  at  inci¬ 
dence  angles  of  0q  =  0°,  30°  and  60°  for  the  half-space  configuration 
and  0q  =  0°  for  the  slab  geometry.  The  normalized  slab  thickness  xq 
varied  from  0.2  to  20.  All  intensity  curves  were  normalized  by  set¬ 
ting  S(  /4tt  equal  to  unity. 

Forward  and  backward  diffuse  intensities  decayed  slower  when  the 

albedo  was  larger  and  incidence  angle  smaller.  The  strongest  diffuse 

intensity  in  the  forward  direction  occurred  near  t  equal  to  unity 

while  in  backward  direction  the  strongest  intensity  occurred  at  the 

air-forest  interface  x  =  0.  The  maximum  forward  diffuse  intensities 

flowed  in  directions  near  0  =  90°  for  small  x.  As  x  and  0  increased 

o 

the  maximum  forward  diffuse  intensities  shifted  away  from  the  verti¬ 
cal.  The  x  dependence  was  much  more  significant  than  the  0Q  depen¬ 
dence.  We  found  that  for  small  x,  incidence  angle  0  did  not  affect 
the  maximum  forward  diffuse  intensity.  At  larger  values  of  x  we 
found  that  the  angular  distribution  of  intensity  evolved  into  a  nar¬ 
row  beam  of  radiation  directed  in  the  0=0°  direction.  For  smaller 
albedos,  beam  narrowing  occurred  at  shorter  penetration  depths  when 
compared  to  larger  values  of  W. 

With  regard  to  the  backward  diffuse  intensity,  a  maximum 

occurred  in  the  0  =  90°  direction  and  decreased  to  a  minimum  value 

in  the  0  =  180°  direction.  For  small  x  and  small  0  the  relative 

o 

backward  diffuse  intensity  was  seen  to  be  essentially  independent  of 
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x •  For  small  qq  the  relative  backward  diffuse  intensity  at  small  x 
was  stronger  than  at  large  x •  However,  at  large  0q  ,  the  reverse 
occurred . 

Power  flux  density  curves  in  the  forward  and  backward  directions 
were  presented  for  various  values  of  albedo  (from  a  low  of  0.5  to  a 
high  of  0.95)  and  for  two  angles  of  incidence  (8  =  0°  and  45°). 
Characteristic  of  such  curves  was  that  for  sufficiently  wide  slabs 
(tq>1)  the  forward  flux  density  F^  quickly  reached  a  maximum  value 
at  x=l)  and  thereafter  exponentially  decayed  as  x  became  larger. 

Larger  albedos  (stronger  scattering)  were  seen  to  yield  slower  atten¬ 
uation  rates  for  F+.  The  backward  flux  density  was  observed  to  reach 
its  maximum  at  x  =  0  as  expected.  Dependence  of  flux  on  incidence 
angles  was  also  examined.  Its  effect  was  shown  to  be  pronounced  for 
smaller  albedos. 

VII.  REFERENCES 

ASt  Abramowitz,  M.  and  Stegun,  I. A.  (eds.)  "Handbook  of  Mathemati¬ 
cal  Functions,"  Dover  Publications,  Inc.,  NY,  1965. 

CZw  Case,  K.  and  Zweifel,  P.,  "Linear  Transport  Theory,"  Addison- 
Wesley,  1969. 

CH  Chandrasekhar,  S.,  "Radiative  Transfer,"  Dover  Publications, 
1960. 

Isl  Ishimaru,  A.,  "Wave  Propagation  and  Scattering  in  Random  Media," 
Vol.  1  (Single  Scattering  and  Transport  Theory),  Academic  Press, 
1978 

Oz  Ozisik,  M.N. ,  "Radiative  Transfer  and  Interactions  with  Con¬ 
duction  and  Convection,"  John  Wiley  and  Sons,  1973. 

VE  Violette,  E.J.,  Espeland,  R.H.,  Schwering,  F.K.,  Whitman,  G.M. , 
"Millimeter  Wave  Propagation  in  Vegetation,"  The  Ninth  DARPA 
Tri-Service  Millimeter-Wave  Conference,  Huntsville,  AL, 

20-22  Oct  1981. 

SJ  Schwering,  F.K.,  Johnson,  R.A.,  Rokkos,  N. ,  Whitman,  G.M., 

"Effects  of  Vegetation  and  Battlefield  Obscurants  on  Point-to- 
Point  Transmission  in  the  Lower  Millimeter  Wave  Region  (30-60 
GHz),"  US  Army  Science  Conference,  West  Point,  NY,  15-18  Jun 
1982. 


25 


ASc  Ziering,  S.  and  Schiff,  D.  ,  "Yvon's  Method  for  Slabs,"  Nuclear 
Science  and  Engineering,  Vol -  3,  pp.  635-647. 

Sny  Snyder,  M.A. ,  "Chebyshev  Methods  in  Numerical  Approximation." 

BFR  Burden,  L.R.,  Faires,  J.D.,  Reynolds,  C.A. ,  "Numerical  Analy¬ 
sis,"  Pr indie,  Wieber  and  Schmidt,  Boston,  1981. 

VIII.  ACKNOWLEDGEMENT 

This  study  was  conducted  jointly  by  the  Center  for  Communi¬ 
cation  Systems  of  the  US  Army  Communications-Electronics  Command  and 
the  New  Jersey  Institute  of  Technology. 

During  the  summers  of  1980  and  1981  Professor  Whitman  worked  under 
the  Laboratory  Research  Cooperative  Program  with  Dr.  F.  Schwering  of 
the  Command's  Millimeter  Wave  and  Secure  Fiber  Team.  Mr.  Li-Wen  Chen 
received  his  MSEE  degree  from  NJIT  in  October  1982  based  on  this  re¬ 
search  and  is  currently  a  professor  at  Koashung  College,  Taiwan. 

The  authors  wish  to  express  their  gratitude  to  Dr.  J.  Robert 
Christian,  Leader  of  the  Millimeter  Wave  Team  and  Dr.  Larry  Dworkin, 
Chief  of  the  Multichannel  Transmission  Division,  for  their  encourage¬ 
ment  and  support. 


26 


Figure  1. 

SLAB  MODEL  OF  A  FOREST  ILLUMINATED  BY  A  PLANE  WAVE. 


Figure  2. 

SPECIFIC  INTENSITY  Ili-A)  FLOWING  THROUGH  AREA  ELEMENT  dA1 


27 


CD 

CD 


® 

CD 

<7» 

<s 

m 

00 

CD 

CD 


CD 

CD 

• 

*0 

CD 

CD 

• 

IO 

CD 

CD 


CD 

CD 

• 

ro 

CD 

CD 

• 

CVJ 

CD 

CD 


CD 

CD 


II 

a 


t 


u-ocaco  fii*-' Luu.  hzhlj2cii*hk> 


oco 


Hi 

o 

<r 

Q. 

CO 

CM 

\ 


OC 

O 

Lu 

Ul 

O 

2E 

H* 

CO 

M 

o 


CO 

I> 


to 

UJ 


LlJ 

CO 


LU 

M 

o 

o 

oc 

<x 

oc 

o 


(0 

O 

3 

Cr> 

■  r~i 

tu 


29 


0.25,  Q0»  30.00 
0.50,  Q0»  30.00° 
0.75,  Q0=»  30.00° 
0.90,  Q0s  30.08P 
0.95,  G0*  30.00° 


U.tt2  & 


0.25,  Q0«  60.00 
0.50,  Q0=  60.00 
0.75,  Q0»  60.00 
0.90,  Q0«  60.00 
0.95,  Q0*  60.00 


© 

© 


© 

© 

© 

© 

© 

© 

© 

t 

© 

© 

© 

© 

© 

© 

© 

© 

© 

• 

• 

• 

0 

• 

• 

» 

• 

© 

© 

© 

© 

© 

© 

© 

© 

<M 

CVI 

ro 

ID 

« 

( 

i 

( 

1 

(JL.lk.ZCl 

Ci  MU.  it 

z  «/> 

Ci© 

31 


TAU — > 

Figure  4.c.  FORWARD  DIFFUSE  INTENSITY  US.  DISTANCE  FOR  1^2  SPACE  WITH  Q=0 
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Figure  5. a.  BACKWARD  DIFFUSE  INTENSITY  US.  DISTANCE 
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Figure  s.b  BACKNARD  DIFFUSE  INTENSITY  US.  DISTANCE  FOR  1/2  SPACE  MITH  0*18^ 
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Figure  5.c.  BACKWARD  DIFFUSE  INTENSITY  US.  DISTANCE  FOR  1/2  SPACE  WITH  Q*180 
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Figure  6.d.  RELATIVE  FRMP  DIFF  INTENSITY  US.  SCATTER  ANGLE  FOR  1/2  SPACE, M 


0.05,  Q0*  30.00 
0.10,  Q0=»  30.00* 
0.20,  Q0*  30.00* 


.  \ 
•.  \ 
%  \ 


o  o  o  o  o 

CD  O 
B  © 

m  m 

CD  CD 

ro  ro 

B  U 
B  B 

©  o 


CD  B 

in  © 

a  a 

©  a-. 


II  H  M  B  II 


— •  <vi  ro  in 


I 

\  / 

v  ... 


B 

m 

B 

in 

B 

B 

IS- 

in 

OJ 

B 

• 

• 

• 

• 

• 

Ol 

B  <X 
B  h- 
•  Ui 
B  X 


B 
B  <r» 
B 
% 

B 


11.20  Njh  S  D-LUQSh- 


RELATIUE  FRWD  DIFF  INTENSITY  MS.  SCA1 


00*09  *00  *0S*0 
00*09  *00  ‘02*0 
00*09  *00  <01 *0 
00*09  *00  <50*0 


VO 

HZH  a, 

u 

3 

CT> 

'H 


40 


<—  THETA  — > 

RELATIVE  FRMD  DIFF  INTENSITY  VS.  SCATTER  ANGLE  FOR  UZ  SPACE, H*0. 75 
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Figure  6.h.  RELATIUE  FRMD  DIFF  INTENSITY  YS.  SCATTER  ANGLE  FOR  1/2  SPACE, W»0.75 
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Figure  6.i.  RELATIUE  FRWD  DIFF  INTENSITY  US.  SCATTER  ANGLE  FOR  1/2  SPACE, H=9. 75 
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igure  6.1-  RELATIVE  FRND  DIFF  INTENSITY  MS.  SCATTER  ANGLE  FOR  UZ  SPACE. H-0. 90 
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Figure’ 7. b.  RELATIVE  FRND  DIFF  INTENSITY  US.  SCATTER  ANGLE  FOR  1/2  SPACE,W=0.2 
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RELATIVE  FRMO  OIFF  INTENSITY  VS.  SCATTER  ANGLE  FOR  1/2  SPACE? M=8. 90 
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Figure  8. a.  RELATIVE  BWD  DIFF  INTENSITY  VS.  SCATTER  ANGLE  FOR  1,'Z  SPACE, W=8 
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Figure  8. C.  RELATIUE  BWD  DIFF  INTENSITY  US.  SCATTER  ANGLE  FOR  1/2  SPACE,U=0.2 


Figure  8. e.  RELATIUE  BWD  DIFF  INTENSITY  US.  SCATTER  ANGLE  FOR  1/2  SPACE, M=0. 50 
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Figure  9. a.  FORWARD  DIFFUSE  INTENSITY  US.  DISTANCE  FOR  SLAB, Q=0°,to=O.2 
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Figure  9. b  FORWARD  DIFFUSE  INTENSITY  US.  DISTANCE  FOR  SLAB, Q=6°» V0  5 
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Figure  9.C.  FORWARD  DIFFUSE  INTENSITY  US.  DISTANCE  FOR  SLAB, Q=0°,t=1  . 
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Figure  9.e  FORWARD  DIFFUSE  INTENSITY  US.  DISTANCE  FOR  SLAB,  G=0°,t  =5.0 
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Figure  10. a  BACKWARD  DIFFUSE  INTENSITY  US.  DISTANCE  FOR  SLAB. Q=18CP,t  =0.2. 
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Figure  10.C.  BACKWARD  DIFFUSE  INTENSITY  US.  DISTANCE  FOR  SLAB, Q=!Stf,T  =1 
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Figure  10. e.  BACKWARD  DIFFUSE  INTENSITY  US.  DISTANCE  FOR  SLAB, Q=18rf*V5  0 
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Figure  11. e  RELATIVE  FWD  DIFF  INTENSITY  VS.  SCATTER  ANGLE  FOR  SLAB,N=0. 
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Figure  ii. g  RELATIVE  FMD  DIFF  INTENSITY  VS.  SCATTER  ANGLE  FOR  SLAB,W=0 
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Figure  12. g.  RELATIUE  BHD  DIFF  INTENSITY  US.  SCATTER  ANGLE  FOR  SLABjQ0=O,W= 
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Figure  I2.h.  RELATIVE  BUD  DIFF  INTENSITY  VS.  SCATTER  ANGLE  FOR  SLAB, Q0=8, W=0. ? 
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Figure  12.  j  RELATIUE  BUD  L»IFF  INTENSITY  US.  SCATTER  ANGLE  FOR  SLAB,  QQ=0°.  W=0. 50 
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12k  RELhTIUE  BUD  DIFF  INTENSITY  US.  SCATTER  ANGLE  FOR  SLAB,Q0=0,N= 
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Figure  13. b.  BACKWARD  DIFFUSE 
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Figure  14. a.  FORWARD  DIFFUSE  FLUX  US.  DISTANCE  FOR  SLAB,T0=  0.25, 
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Figure  14. b.  BACKWARD  DIFFUSE  FLUX  US.  DISTANCE  FOR  SLAB,T0=  0.2 


Figure  15. a.  FORWARD  DIFFUSE  FLUX  US.  DISTANCE  FOR  SLAB,TO=  1.80,S= 


00*0  =00 


CD 

© 

in 


CD  ©  © 
©  ©  © 
•  m  « 

©  ©  © 


II  II  II 
Q  (J)  QD 

G  G  G 

©  in  ©  in 
in  cri  <x» 

©  ©  ©  © 


ii  ii  ii  ii 

2  3  3  3 


— «  cu  ro  rf 


©. 

© 

© 

© 

©. 

© 

© 

©* 

©* 

©* 

© 

©’ 

© 

ro 

'vO 

<T> 

CM 

in 

i 

i 

1 

w—4 

© 

© 

•*r 


© 

© 

ro 


© 

© 

CM 


© 

© 


© 

© 

© 


t 

3 

<r 


© 

© 

in 


CO  ^  3  Cl  OMU.U.  U - OX  •—  2  Cl© 


124 


Figure  i6.b.  BACKWARD  DIFFUSE  FLUX  US.  DISTANCE  FOR  SLAB,T6= 
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Figure  17. a.  FORWARD  DIFFUSE  FL 
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Figure  17. b.  BACKMhRD  DIFFUSE  FLUX  US.  DISTANCE  FOR  SLAB,  T0=2G.  O0, 
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Figure  19. a.  FORWARD  FLUX  US.  DISTANCE  FOR  SLhE:,TQ= 
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Figure  20. a.  FORWARD  FLUX  US.  DISTANCE  FOR  SLAB,T0=20. 06, 
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Figure  20. b.  BACKWARD  DIFFUSE  FL 
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Figure  2i.a.  FORWARD  DIFFUSE  FLUX  US.  DISTANCE, S#= 
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Figure  21.c.  BACKWARD  DIFFUSE 
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Figure  2i.e.  BACKWARD  DIFFUSE  FLUX  US.  DISTANCE,  SLAB, 


TABLE  1 . 

EIGENVALUES  FOR  W  =  0.50,  N  =  15 


Eigenvalues 
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TABLE  2. 

EIGENVALUES  FOR  W=0.75,  N  =  15 


Eigenvalues 


N= 15, W=0. 75 


0.0024151615 


0.0217681920 


0.0600775240 


0.1162469600 


0.1882330600 


0.4694660100 


0.5723481100 


0.6727326500 


0.7664521500 


0.8490449300 


0.9171296600 


0.9663557400 


0.9949324100 


1.2894637000 


Asymptotic 

Eigenvalue 


«•  M 

- -  —  ~ 

0.2734405900 

' 

M  •• 

0.3684320600 

1.2894634000 
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TABLE  3. 

EIGENVALUES  FOR  W  =  0.95,  N  =  15 


Eigenvalues 


Asymptotic 

Eigenvalue 


0.0024176029 


0.0218328890 


0.0603619480 


0.1170050200 


0.1897284200 


0.2759308300 


0.3720490100 
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2.6351487000 


2.6351488000 


